Variance quantitative trait loci mapping in the MoBa cohort: detecting interaction effects involved in childhood behavioural traits

Supplement 1: Analysis of behaviour

Table 1

Sample age and sex distributions of MoBa children (with genetic data and >75% of outcomes)

Df <- fread("N:/durable/projects/crayner/PipePsy/Data/ChPsyFsc.tsv") %>% select(PrID,starts_with("Ch"))
names(Df) <- str_replace_all(names(Df),"Ch","")
names(Df) <- str_replace_all(names(Df),"Q","_Q")
names(Df) <- str_replace_all(names(Df),"Q04","6mo")
names(Df) <- str_replace_all(names(Df),"Q05","18mo")
names(Df) <- str_replace_all(names(Df),"Q06","3yr")
names(Df) <- str_replace_all(names(Df),"Q07","5yr")
names(Df) <- str_replace_all(names(Df),"Q09","8yr")
names(Df) <- str_replace_all(names(Df),"_Q11C","C_14yr")
names(Df) <- str_replace_all(names(Df),"_Q11M","M_14yr")
Df <- Df %>% select(PrID,matches("6mo"),matches("18mo"),matches("3yr"),matches("5yr"),matches("8yr"),matches("14yr"))

excl <- c(
  "IcqFull_6mo","CbclFull_18mo","EasFull_18mo","CbclFull_3yr","CfqFull_3yr",
  "EasFull_3yr","ItseaFull_3yr","ScqFull_3yr","SdqFull_3yr","CbclFull_5yr",
  "CccFull_5yr","EasFull_5yr","SprakFull_5yr","CccFull_8yr","NhpicFull_8yr",
  "RsdbdFull_8yr","ScqFull_8yr","SprakFull_8yr","PeqFullC_14yr","SbFullC_14yr",
  "SclFullC_14yr","SdqFullC_14yr","SdqFullM_14yr","SppaFullC_14yr",
  "CbclAdhd_18mo","CbclAdhd_3yr","CbclAdhd_5yr"
)

Df <- Df %>% select(!c(all_of(excl)))

Cov <- 
  # data.table::fread("N:/durable/projects/crayner/Data/Participation/MoBaAgeVarsLong.csv")%>% 
  data.table::fread("N:/durable/projects/crayner/Data/Participation/MoBaAgeVarsLongLinked.csv") %>% 
  select(MoLNR,FaLNR,PrID,ChIID,ChN,ChAge,MoAge,FaAge,Wave=Q) %>% # ChSex=ChSEX,
  filter(ChAge>0&Wave!="QF2"&Wave!=""&Wave!=" ") %>% 
  distinct() %>% # filter(ChAge>0) %>% 
  mutate(PrID=paste0(PrID,"_",ChN), Wave=droplevels(as.factor(Wave))) %>% 
  filter(PrID %in% Df$PrID) %>% 
  mutate(FamID=as.numeric(as.factor(paste0(MoLNR,"_",FaLNR)))) %>% 
  select(PrID,FamID,ChIID,Wave,matches("Age")) %>% #,ChSex
  mutate(Wave=factor(Wave,levels=c("Q04","Q05","Q06","Q07","Q09","Q11C","Q11M"),
                     labels=c("6mo","18mo","3yr","5yr","8yr","C14yr","M14yr"))) %>% 
  filter(ChIID!=""&ChIID!=" "&!is.na(ChIID))

AgeDiff <- Cov %>% filter(Wave=="6mo") %>% mutate(ChMo=MoAge-ChAge,ChFa=FaAge-ChAge) %>%  select(PrID,ChMo,ChFa) 

Cov2 <- data.table::fread("N:/durable/projects/crayner/Data/Covariates/ChAge14yrLong.tsv") %>% 
  inner_join(AgeDiff) %>% 
  mutate(MoAge=ChMo+ChAge,FaAge=ChFa+ChAge)%>%select(-ChMo,-ChFa) 

Cov3 <- data.table::fread("N:/durable/projects/crayner/Data/Covariates/MbrnCovariates113536.tsv")%>% 
  mutate(PrID=paste0(PrID,"_",BaN)) %>% 
  filter(PrID %in% Df$PrID) %>% 
  mutate(ChSex=factor(ChildSGender,levels=c("Male","Female"),labels=c(1,2))) %>% 
  select(PrID,ChSex) 

Cov <- Cov %>% bind_rows(Cov2) %>% inner_join(Cov3)

tmpDf_1 <- 
  Df %>% 
    select(PrID,matches("6mo")) %>% 
    mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.1) %>% #na.omit() %>% 
    select(PrID) %>% inner_join(Cov %>% filter(Wave=="6mo")) %>% 
    summarize( N=n(), 
             `Males`=sum(ChSex==1,na.rm=T),
             `Males %`=sum(ChSex==1,na.rm=T)/n(),
             `Females`=sum(ChSex==2,na.rm=T), 
             `Females %`=sum(ChSex==2,na.rm=T)/n(), 
             `Child Age (mean)`=round(mean(ChAge),2), 
             `Mothers Age (mean)`=round(mean(MoAge),2), 
             `Fathers Age (mean)`=round(mean(FaAge),2),
             `Child Age (sd)`=round(sd(ChAge),2), 
             `Mothers Age (sd)`=round(sd(MoAge),2), 
             `Fathers Age (sd)`=round(sd(FaAge),2)) %>%              
    select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t()%>% as.data.frame()%>% rownames_to_column("labs")

tmpDf_2 <- 
  Df %>% 
    select(PrID,matches("18mo")) %>% 
    mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.35) %>% #na.omit() %>% 
    select(PrID) %>% inner_join(Cov %>% filter(Wave=="18mo")) %>% 
    summarize( N=n(), 
             `Males`=sum(ChSex==1,na.rm=T),
             `Males %`=sum(ChSex==1,na.rm=T)/n(),
             `Females`=sum(ChSex==2,na.rm=T), 
             `Females %`=sum(ChSex==2,na.rm=T)/n(), 
             `Child Age (mean)`=round(mean(ChAge),2), 
             `Mothers Age (mean)`=round(mean(MoAge),2), 
             `Fathers Age (mean)`=round(mean(FaAge),2),
             `Child Age (sd)`=round(sd(ChAge),2), 
             `Mothers Age (sd)`=round(sd(MoAge),2), 
             `Fathers Age (sd)`=round(sd(FaAge),2)) %>%              
    select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t()%>% as.data.frame()%>% rownames_to_column("labs")

tmpDf_3 <- 
  Df %>% 
    select(PrID,matches("3yr")) %>% 
    mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.1) %>% #na.omit() %>% 
    select(PrID) %>% inner_join(Cov %>%  filter(Wave=="3yr")) %>% 
    summarize( N=n(), 
             `Males`=sum(ChSex==1,na.rm=T),
             `Males %`=sum(ChSex==1,na.rm=T)/n(),
             `Females`=sum(ChSex==2,na.rm=T), 
             `Females %`=sum(ChSex==2,na.rm=T)/n(), 
             `Child Age (mean)`=round(mean(ChAge),2), 
             `Mothers Age (mean)`=round(mean(MoAge),2), 
             `Fathers Age (mean)`=round(mean(FaAge),2),
             `Child Age (sd)`=round(sd(ChAge),2), 
             `Mothers Age (sd)`=round(sd(MoAge),2), 
             `Fathers Age (sd)`=round(sd(FaAge),2)) %>%              
    select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t()%>% as.data.frame()%>% rownames_to_column("labs")

tmpDf_4 <- 
  Df %>% 
    select(PrID,matches("5yr")) %>% 
    mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.1) %>% #na.omit() %>% 
    select(PrID) %>% inner_join(Cov %>%  filter(Wave=="5yr")) %>% 
    summarize( N=n(),
             `Males`=sum(ChSex==1,na.rm=T),
             `Males %`=sum(ChSex==1,na.rm=T)/n(),
             `Females`=sum(ChSex==2,na.rm=T), 
             `Females %`=sum(ChSex==2,na.rm=T)/n(), 
             `Child Age (mean)`=round(mean(ChAge),2), 
             `Mothers Age (mean)`=round(mean(MoAge),2), 
             `Fathers Age (mean)`=round(mean(FaAge),2),
             `Child Age (sd)`=round(sd(ChAge),2), 
             `Mothers Age (sd)`=round(sd(MoAge),2), 
             `Fathers Age (sd)`=round(sd(FaAge),2)) %>%              
    select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t()%>% as.data.frame()%>% rownames_to_column("labs")

tmpDf_5 <- 
  Df %>% 
    select(PrID,matches("8yr")) %>% 
    mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.1) %>% #na.omit() %>% 
    select(PrID) %>% inner_join(Cov %>%  filter(Wave=="8yr")) %>% 
    summarize( N=n(), 
               `Males`=sum(ChSex==1,na.rm=T),
               `Males %`=sum(ChSex==1,na.rm=T)/n(),
               `Females`=sum(ChSex==2,na.rm=T), 
               `Females %`=sum(ChSex==2,na.rm=T)/n(), 
               `Child Age (mean)`=round(mean(ChAge),2), 
               `Mothers Age (mean)`=round(mean(MoAge),2), 
               `Fathers Age (mean)`=round(mean(FaAge),2),
               `Child Age (sd)`=round(sd(ChAge),2), 
               `Mothers Age (sd)`=round(sd(MoAge),2), 
               `Fathers Age (sd)`=round(sd(FaAge),2)) %>%              
    select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t()%>% as.data.frame()%>% rownames_to_column("labs")

tmpDf_6 <- 
  Df %>% 
  select(PrID,matches("C_14yr")) %>% 
  mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.5) %>% #na.omit() %>% 
  select(PrID) %>% inner_join(Cov %>% filter(Wave=="C14yr")) %>%  #) %>% distinct(PrID,.keep_all=T)
  summarize( N=n(),
             `Males`=sum(ChSex==1,na.rm=T),
             `Males %`=sum(ChSex==1,na.rm=T)/n(),
             `Females`=sum(ChSex==2,na.rm=T), 
             `Females %`=sum(ChSex==2,na.rm=T)/n(), 
             `Child Age (mean)`=round(mean(ChAge),2), 
             `Mothers Age (mean)`=round(mean(MoAge),2), 
             `Fathers Age (mean)`=round(mean(FaAge),2),
             `Child Age (sd)`=round(sd(ChAge),2), 
             `Mothers Age (sd)`=round(sd(MoAge),2), 
             `Fathers Age (sd)`=round(sd(FaAge),2)) %>%              
  select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t() %>% as.data.frame() %>% rownames_to_column("labs")

tmpDf_7 <- 
  Df %>% 
  select(PrID,matches("M_14yr")) %>% 
  mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% na.omit() %>% #na.omit() %>% 
  select(PrID) %>% inner_join(Cov %>% filter(Wave=="M14yr")) %>%  #) %>% distinct(PrID,.keep_all=T)
  summarize( N=n(),
             `Males`=sum(ChSex==1,na.rm=T),
             `Males %`=sum(ChSex==1,na.rm=T)/n(),
             `Females`=sum(ChSex==2,na.rm=T), 
             `Females %`=sum(ChSex==2,na.rm=T)/n(), 
             `Child Age (mean)`=round(mean(ChAge),2), 
             `Mothers Age (mean)`=round(mean(MoAge),2), 
             `Fathers Age (mean)`=round(mean(FaAge),2),
             `Child Age (sd)`=round(sd(ChAge),2), 
             `Mothers Age (sd)`=round(sd(MoAge),2), 
             `Fathers Age (sd)`=round(sd(FaAge),2)) %>%              
  select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t() %>% as.data.frame() %>% rownames_to_column("labs")

tmpDf <- plyr::join_all(list(tmpDf_1,tmpDf_2,tmpDf_3,tmpDf_4,tmpDf_5,tmpDf_6,tmpDf_7),by="labs") %>% column_to_rownames("labs")
names(tmpDf)<- c("6mo","18mo","3yr","5yr","8yr","C14yr","M14yr")

datatable(
  tmpDf %>% mutate(across(everything(), function(x) round(x,2))),
  extensions='Buttons',
  options=list(
    dom='tiB',buttons=c('csv'),
    searching=F,paging=F,ordering=F,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

Time-points named to reflect the average age of the children: 6mo, 6 months old; 18mo, 18 months old; 3yr, 3 years old; 5yr, 5 years old; 8yr, 8 years old; 14yr, 14 years old. Outcomes measured at all time-points are mother-rated, apart from at the 14yr time-point, where there are also child ratings, indicated by c.

S.Table 1

Psychometric outcomes measured in MoBa children

Df <- fread("N:/durable/projects/crayner/PipePsy/Data/ChPsyFsc.tsv") %>% select(PrID,starts_with("Ch"))
names(Df) <- str_replace_all(names(Df),"Ch","")
names(Df) <- str_replace_all(names(Df),"Q","_Q")
names(Df) <- str_replace_all(names(Df),"Q04","6mo")
names(Df) <- str_replace_all(names(Df),"Q05","18mo")
names(Df) <- str_replace_all(names(Df),"Q06","3yr")
names(Df) <- str_replace_all(names(Df),"Q07","5yr")
names(Df) <- str_replace_all(names(Df),"Q09","8yr")
names(Df) <- str_replace_all(names(Df),"_Q11C","C_14yr")
names(Df) <- str_replace_all(names(Df),"_Q11M","M_14yr")
Df <- Df %>% select(PrID,matches("6mo"),matches("18mo"),matches("3yr"),matches("5yr"),matches("8yr"),matches("14yr"))

excl <- c(
  "IcqFull_6mo","CbclFull_18mo","EasFull_18mo","CbclFull_3yr","CfqFull_3yr",
  "EasFull_3yr","ItseaFull_3yr","ScqFull_3yr","SdqFull_3yr","CbclFull_5yr",
  "CccFull_5yr","EasFull_5yr","SprakFull_5yr","CccFull_8yr","NhpicFull_8yr",
  "RsdbdFull_8yr","ScqFull_8yr","SprakFull_8yr","PeqFullC_14yr","SbFullC_14yr",
  "SclFullC_14yr","SdqFullC_14yr","SdqFullM_14yr","SppaFullC_14yr",
  "CbclAdhd_18mo","CbclAdhd_3yr","CbclAdhd_5yr"
)

Df <- Df %>% select(!c(all_of(excl)))

# excl <- names(Df)[grepl(x=names(Df),"Full")]
# fullScales <- c(
#   "CastFull_5yr","CprsrsFull_5yr","SlasFull_5yr","CebqFull_8yr","ScaredFull_8yr",
#   "SmfqFull_8yr","SpinFull_8yr","CapeFullC_14yr","SmfqFullC_14yr","SwlsFullC_14yr",
#   "YptiFullC_14yr"
# )
# fullScales <- c("Cast","Cprsrs","Slas","Cebq","Scared","Smfq","Spin","Cape","Smfq","Swls","Ypti")
# test <- names(Df)[grepl(x=names(Df),paste0(fullScales,collapse="|"))]

Tab1 <- lapply(Df, function(x) sum(!is.na(x)))
Tab1 <- do.call(Tab1,what="bind_rows") 
Tab1 <- t(Tab1) %>% as.data.frame()
Tab1 <- Tab1 %>% rownames_to_column("Var")
Tab1 <- Tab1[-1,]
Tab1 <- Tab1 %>% mutate(Var=str_replace_all(Var,"C_14yr", "_14yr Child"))
Tab1 <- Tab1 %>% mutate(Var=str_replace_all(Var,"M_14yr", "_14yr Mother"))
Tab1 <- Tab1 %>% mutate(Var=str_replace_all(Var,"([a-z])([A-Z])", "\\1 \\2"))
Tab1 <- Tab1 %>% mutate(Var=str_replace_all(Var,"_", " "))
Tab1 <- Tab1 %>% separate(Var, c("Scale","Subscale","Age","Rater"))
Tab1 <- Tab1 %>% mutate(Rater=ifelse(is.na(Rater),"Mother",Rater))
Tab1 <- Tab1 %>% mutate(Scale=str_to_upper(Scale))

datatable(
  Tab1 %>% select(Scale,Subscale,Age,Rater,N=V1),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

ICQ6, Infant Characteristics Questionnaire 6 Months Form; ASQ, Ages and Stages Questionnaires; EAS, Emotionality, Activity and Shyness Temperament Questionnaire; CBCL, Child Behaviour CheckList; AUTQ, Early Screening of Autistic Traits Questionnaire-Short (+Modified Checklist for Autism in Toddlers); ITSEA, Infant-Toddler Social and Emotional Assessment; CFQ, The Child Feeding Questionnaire; SDQ, Strength and Difficulties Questionnaire-Prosocial Subscale; SCQ, Social Communication Questionnaire; PPBS, Preschool Play Behaviour Scale; CPRSRS, Conners Parent Rating Scale - Revised, Short Form; CCC, Children’s Communication Checklist-2 Short: Coherence I; SPRAK, Statements about Language-Related Difficulties: Semantics; SLAS, Speech and Language Assessment Scale; CAST, Childhood Asperger Syndrome Test; RSDBD, Rating Scale for Disruptive Behaviour Disorders; CEBQ, Children’s Eating Behaviour Questionnaire; SCARED, Screen for Child Anxiety Related Disorders; SMFQ, Short Mood and Feelings Questionnaire; SPIN, Mini Social Phobia Inventory; NHPIC30, Short Norwegian Hierarchical Personality Inventory for Children; PEQ, Parental Environment Questionnaire; YPTI, Youth Psychopathic Traits; CAPE15, Community Assessment of Psychic Experiences; SCL, The (Hopkins) Symptoms Checklist; DES, Differential Emotional Scale; RSES, Rosenberg Self-Esteem Scale; SWLS, Satisfaction with life scale; IPIP, International Personality Item Pool; SFS, School Functioning Scale; SPPA, Self-Perception Profile for Adolescents; SB, School belonging

S.Table 2

Item response theory model fit statistics for all psychometric scales

ModStat <- fread("Tables_Phe/xIRT_ModelStats.tsv") %>%
  filter(grepl("Ch",Scale)) %>% 
  mutate(Trait=paste0(Scale,Subscale)) %>% 
  select(Trait,Age=Wave,everything()) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  mutate(Age=str_replace_all(Age,"Q04","6mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q05","18mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q06","3yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q07","5yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q09","8yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>% 
  mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>% 
  select(!matches("Scale"),!matches("Subscale")) %>% 
  filter(!Trait %in% excl) %>% 
  select(-Scale,-p,-RMSEA_5,-RMSEA_95,-LL)


datatable(
  ModStat,
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=T,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,scrollY=T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

1 See table one for description of each trait. 2 Time-points named to reflect the average age of the children: 6mo, 6 months old; 18mo, 18 months old; 3yr, 3 years old; 5yr, 5 years old; 8yr, 8 years old; 14yr, 14 years old. Outcomes measured at all time-points are mother-rated, apart from at the 14yr time-point, where there are also child ratings, indicated by c. 3 The item response theory model which had the best fit to the data: graded, graded response model (ref); nominal, nominal response model; 3-PL, three parameter logistic model (ref); Rasch model (ref). 4AIC, Akaike Information Criterion. 5 BIC, Bayesian Information Criterion. 6 Model fit assessed using M2 index (2 variant), which is specifically designed to assess the fit of item response models for ordinal data - when polytomous response models do not have sufficient degrees of freedom to compute M2. This function also computes associated fit indices that are based on fitting the null model. 7 The degrees of freedom for the M2 test. 8 RMSEA, the M2-based root mean square error of approximation is the primary fit index - can be used to suggest that the data fit the model, with suggested cutoff values of RMSEA <= .06. 9 SRMSR, standardised root mean square residual and 10 CFI, comparative fit index is used to assess adequacy of model fit. The suggested cutoff values of SRMSR <= .08. 12 Reliability, a measure of the average inter-item correlation, assessed using Cronbach’s alpha.

S.Figure 2

Phenotypic correlations (a) across trait levels (left) and trait variability (right), with associated scree plots (b) and factor loadings following EFA (c).

A

Factor scores

# ESTIMATING MATRICES (HETCOR)
CorMat   <- cormat(Df %>% select(-PrID))
# Make interactive correlation table
CorTably <- cortab(CorMat)

# Plot interactive correlation matrix
CorTab   <- cortably(CorMat)
CorTab$z <- abs(CorTab$r/CorTab$se)

Cor <- 
 CorTab %>%
 select(Var1,Var2,r) %>%
 filter(Var1!=Var2) %>%
 pivot_wider(names_from = Var2, values_from = r) %>%
 column_to_rownames("Var1") 

Err  <- 
 CorTab %>%
 select(Var1,Var2,se) %>%
 filter(Var1!=Var2) %>%
 pivot_wider(names_from = Var2, values_from = se) %>%
 column_to_rownames("Var1")

Zsc  <- 
 CorTab %>%
 select(Var1,Var2,z) %>%
 filter(Var1!=Var2) %>%
 pivot_wider(names_from = Var2, values_from = z) %>%
 column_to_rownames("Var1")

Zsc  <- log10(Zsc)

ChCorplotlyFull <- heatmaply::heatmaply_cor(
  as.data.frame(Cor),
  node_type = "scatter",
  point_size_mat = Zsc,
  point_size_name = "log10|Z|",
  label_names = c("x", "y", "r"),
  showticklabels = c(T,F),
  fontsize_col=7,
  grid_size = 0.001,
  dendrogram = "both",
  subplot_heights=c(0.04, 0.96),
  subplot_widths=c(0.94,0.06),
  hide_colorbar = T
  )

ChCorplotlyFull$height <- nrow(Cor)*13
ChCorplotlyFull$width  <- nrow(Cor)*13
# ChCorplotlyFull$x$data[[2]]$marker$size <- (ChCorplotlyFull$x$data[[2]]$marker$size/3)
# ChCorplotlyFull$x$data[[3]]$marker$size <- (ChCorplotlyFull$x$data[[3]]$marker$size/3)
ChCorplotlyFull$x$data[[4]]$marker$size <- (ChCorplotlyFull$x$data[[4]]$marker$size/2.5)
# ChCorplotlyFull$x$data[[5]]$marker$size <- (ChCorplotlyFull$x$data[[5]]$marker$size/3)

ChCorplotlyFull

Quantile integrated rank scores

vDf <- fread("N:/durable/projects/crayner/PipePsy/Data/FamPsyQrs.tsv") %>% select(PrID,starts_with("Ch"))
names(vDf) <- str_replace_all(names(vDf),"Ch","")
names(vDf) <- str_replace_all(names(vDf),"Q","_Q")
names(vDf) <- str_replace_all(names(vDf),"Q04","6mo")
names(vDf) <- str_replace_all(names(vDf),"Q05","18mo")
names(vDf) <- str_replace_all(names(vDf),"Q06","3yr")
names(vDf) <- str_replace_all(names(vDf),"Q07","5yr")
names(vDf) <- str_replace_all(names(vDf),"Q09","8yr")
names(vDf) <- str_replace_all(names(vDf),"_Q11C","C_14yr")
names(vDf) <- str_replace_all(names(vDf),"_Q11M","M_14yr")
vDf <- vDf %>% select(PrID,matches("6mo"),matches("18mo"),matches("3yr"),matches("5yr"),matches("8yr"),matches("14yr"))

vDf <- vDf %>% select(!c(excl))

# ESTIMATING MATRICES (HETCOR)
# vCorMat <- cormat(vCorvDf)
vCorMat <- cormat(vDf %>% select(-PrID))
# Make interactive correlation table
vCorTably <- cortab(vCorMat)
# Plot interactive correlation matrix
vCorTab <- cortably(vCorMat)

vCorTab$z <- abs(vCorTab$r/vCorTab$se)


vCor <- 
 vCorTab %>%
 select(Var1,Var2,r) %>%
 filter(Var1!=Var2) %>%
 pivot_wider(names_from = Var2, values_from = r) %>%
 column_to_rownames("Var1") 

vCor[is.na(vCor)] <- 0

Err  <- 
 vCorTab %>%
 select(Var1,Var2,se) %>%
 filter(Var1!=Var2) %>%
 pivot_wider(names_from = Var2, values_from = se) %>%
 column_to_rownames("Var1")

Err[is.na(Err)] <- 0

Zsc  <- 
 vCorTab %>%
 select(Var1,Var2,z) %>%
 filter(Var1!=Var2) %>%
 pivot_wider(names_from = Var2, values_from = z) %>%
 column_to_rownames("Var1")

Zsc  <- log10(Zsc)

#TO DO: ADD SE LABELS
ChvCorplotlyFull <- heatmaply::heatmaply_cor(
  as.data.frame(vCor),
  node_type = "scatter",
  point_size_mat = Zsc,
  point_size_name = "log10|Z|",
  # point_size_mat = logP,
  # point_size_name = "-log10P/100",
  label_names = c("x", "y", "r"),
  showticklabels = c(T,F),
  fontsize_col=7,
  grid_size = 0.001,
  dendrogram = "both",
  subplot_heights=c(0.06, 0.94),
  subplot_widths=c(0.94,0.06),
  hide_colorbar = T
  )

ChvCorplotlyFull$height <- nrow(vCor)*13
ChvCorplotlyFull$width  <- nrow(vCor)*13
# ChvCorplotlyFull$x$data[[2]]$marker$size <- (ChvCorplotlyFull$x$data[[2]]$marker$size/3)
# ChvCorplotlyFull$x$data[[3]]$marker$size <- (ChvCorplotlyFull$x$data[[3]]$marker$size/3)
ChvCorplotlyFull$x$data[[4]]$marker$size <- (ChvCorplotlyFull$x$data[[4]]$marker$size/2.5)
# ChvCorplotlyFull$x$data[[5]]$marker$size <- (ChvCorplotlyFull$x$data[[5]]$marker$size/3)

ChvCorplotlyFull

B

Factor scores

screePlot  <- scree_plot(data=Df%>%select(-PrID), method = "ml", n.iter=200)+ scale_x_continuous(breaks=1:25,limits=c(1,25))
screePlot1 <- screePlot + scale_x_continuous(breaks=1:25,limits=c(1,25))
screePlot1

Quantile integrated rank scores

screePlot  <- scree_plot(data=vDf%>%select(-PrID), method = "ml", n.iter=200)+ scale_x_continuous(breaks=1:25,limits=c(1,25))
screePlot2 <- screePlot + scale_x_continuous(breaks=1:25,limits=c(1,25))
screePlot2

C

Factor scores

Parallel

PLOT1

Kaiser

PLOT2

Quantile integrated rank scores

Parallel

PLOT3

Kaiser

PLOT4

S.Table 3

Phenotypic correlations

Factor scores

datatable(
  CorTably %>% select(Var1,Var2,r,se,p,n),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=T,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,scrollY=T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

Quantile integrated rank scores

datatable(
  vCorTably %>% select(Var1,Var2,r,se,p),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

S.Figure 1

Plots of SNP effects on phenotype mean and variance.

Violin plots display the distribution of a phenotype by minor allele count, at genotypes with different effects on the mean and variance of the phenotype. (A-D) The red dotted lines represent the regression of phenotype on genotype. Blue lines represent regression of phenotype on genotype, restricted to individuals in the upper and lower quintiles of the phenotype distribution - this helps to visualise the heteroscedasticity across genotype levels. All genotypes were simulated with minor allele frequency equal to 0.4 and environments had a normal distribution (n=1000). All phenotypes were simulated as a function of genetic and/or environmental effects with random normal error distribution. (A) For the null model the phenotype was simulated without any effect from genotype or environment. (B) the mQTL model was simulated with genetic effect 𝛽=.2. (C) the vQTL model was simulated with genetic effect 𝛽=0, an environmental effect 𝛽=.5 and a GxE interaction effect 𝛽=.4. (D) the mvQTL model was simulated with genetic effect 𝛽=.2, an environmental effect 𝛽=.5 and a GxE interaction effect 𝛽=.4. (E-F) Points are coloured based on a second variable. Lines represent regression of phenotype on genotype stratified by (E) tertiles of environmental exposure and (F) a genotype at a second locus. This illustrates how the interaction between a genotype and a second factor (an environmental exposure or a second genetic locus) contributes to increased variance across the genotype levels. (E) the GxE model was simulated with genetic effect 𝛽=.2, an environmental effect 𝛽=.5 and a GxE interaction effect 𝛽=.4. (F) the GxG model was simulated with genotype 1 effect 𝛽=.2, genotype 2 effect 𝛽=.2 and a GxG interaction effect 𝛽=.2.

Figure 1

Quantile-quantile (left) and P-P plots (right) for vGWA and mGWA p-values

Quantile-quantile plots (left) display the expected and observed p-values from GWA (+) and vGWA (x). Scatter plots (right) show the relationship between p-values for SNPs estimated from vGWA on the y-axis and GWA on the x-axis. The relationship between p-values for ASQ motor skills is strongly positively correlated (i.e. the same SNPs are identified as having effects), whereas for ITSEA emotional dysregulation SNPs identified as having stronger variance effects did not necessarily have mean effects.

S.Table 4

Genome-wide significant SNPs identified in mGWA & vGWA analyses (independent loci r²<.1)

SnpTab <- fread("Tables_Snps/ChPsy/ClumpedTopSnps.tsv")

SnpTab$af <- round(SnpTab$af,2)
SnpTab$mP <- formatC(SnpTab$mP, format = "e", digits = 2)
SnpTab$vP <- formatC(SnpTab$vP, format = "e", digits = 2)

DT::datatable(
  SnpTab,
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

CorP <- cor.test(as.numeric(SnpTab$mP),as.numeric(SnpTab$vP))

¹ See table one for description of each trait. ² Time-points named to reflect the average age of the children: 6mo, 6 months old; 18mo, 18 months old; 3yr, 3 years old; 5yr, 5 years old; 8yr, 8 years old; 14yr, 14 years old. Outcomes measured at all time-points are mother-rated, apart from at the 14yr time-point, where there are also child ratings, indicated by c.

Correlation between p-values for mGWA and vGWA: Pearson’s r = -0.36; p = 7.16e-01

Figure 2

PheWas associations across all independent (r<.1) MoBa vGWA (top) and GWA (bottom) associations (p<1e-5)

vPheWasTab <- fread("Tables_Snps/ChPsy/vPheWasTab.tsv")
mPheWasTab <- fread("Tables_Snps/ChPsy/mPheWasTab.tsv")
PheWasTab3 <- vPheWasTab %>% bind_rows(mPheWasTab)

PheWasLab <- 
  PheWasTab3 %>% 
  arrange(PheWasP) %>% distinct(PheWasTrait,MoBaTrait,.keep_all=T) %>% 
  select(PheWasTrait,MoBaTrait,PheWasP) %>%
  mutate(Label=str_replace_all(PheWasTrait, "\\s*\\([^\\)]+\\)", "")) %>% 
  mutate(Label=str_replace_all(Label, "[[:punct:]]", "")) %>% 
  mutate(Label=str_to_title(Label)) %>% 
  mutate(Label=str_replace_all(Label, " ", "")) %>% 
  mutate(Label=str_replace_all(Label, "CellCount", "Count")) %>% 
  distinct(Label,MoBaTrait,.keep_all=T) %>% 
  inner_join(PheWasTab3,c("MoBaTrait","PheWasTrait","PheWasP")) %>% 
  select(Label,PheWasTrait,MoBaTrait,vP,mP,PheWasP,rsid) %>% distinct() %>% 
  mutate(MoBaTrait=str_replace_all(MoBaTrait, "C_", " ")) %>% 
  mutate(MoBaTrait=str_replace_all(MoBaTrait, "Full", "")) %>% 
  # mutate(MoBaTrait=str_replace_all(MoBaTrait, "([a-z])([A-Z])", "\\1 \\2")) %>% 
  mutate(MoBaTrait=str_replace_all(MoBaTrait, "([a-z])([0-9])", "\\1 \\2"))%>% 
  separate(MoBaTrait,c("MoBaTrait","Age")," ")%>% 
  mutate(Age = factor(Age, levels=c("6mo","18mo","3yr","5yr","8yr","14yr"))) #%>% 
  # mutate(MoBaTrait=paste0(Scale,Trait))

PheWasLab <- 
  bind_rows(
  PheWasLab %>% filter(vP<5e-8) %>% rename("GwaP"=vP) %>% mutate(Effect="vQTL"),
  PheWasLab %>% filter(mP<5e-8) %>% rename("GwaP"=mP) %>% mutate(Effect="mQTL")
  ) %>% 
  mutate(GwaP=round(-log10(GwaP),1),PheWasP=round(-log10(PheWasP),1))


PheWasTab <- PheWasTab3 %>% select(rsid,PheWasTrait,MoBaTrait,PheWasP,mP,vP)

PheWasTab$PheWasP <- formatC(PheWasTab$PheWasP, format = "e", digits = 2)
PheWasTab$mP <- formatC(PheWasTab$mP, format = "e", digits = 2)
PheWasTab$vP <- formatC(PheWasTab$vP, format = "e", digits = 2)

PheWasPlot <- 
  ggplot(
    PheWasLab,aes(y=PheWasP,x=GwaP,shape=MoBaTrait,colour=MoBaTrait,label=PheWasTrait))+ 
    geom_point()+theme_minimal()+
    facet_grid(cols=vars(Age),rows=vars(Effect))+
    scale_colour_viridis_d(begin=.2,end=.8)+
    scale_shape_manual(values=c(0:14))+
    theme(legend.position="bottom")

library(plotly)

ggplotly(PheWasPlot, tooltip = c("MoBaTrait", "GwaP", "PheWasTrait", "PheWasP"))

# PheWasPlotly <- plotly::ggplotly(PheWasPlot)
# saveWidget(PheWasPlotly, "PsyPheWasPlotly.html", selfcontained = T, libdir = "lib")

Scatter plots display the -log₁₀ GWA p-values from publicly available datasets (PheWasP; y-axis) and -log10 GWA p-values from MoBa analyses (x-axis). The plot is stratified by Age of MoBa participants, with age increasing from left to right (x-strip) and by the type of analysis (y-strip), with mGWA in the top panel and vGWA in the bottom panel. MoBa traits are labelled by the shape and the colour of the point. PheWas traits are labelled in hover over boxes (interactive version) and can be viewed in S.Table 5.

S.Table 5

Results from a phenome wide association scan of publicly available datasets

DT::datatable(
  PheWasTab,
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

¹ See table one for description of each trait. ² Time-points named to reflect the average age of the children: 6mo, 6 months old; 18mo, 18 months old; 3yr, 3 years old; 5yr, 5 years old; 8yr, 8 years old; 14yr, 14 years old. Outcomes measured at all time-points are mother-rated, apart from at the 14yr time-point, where there are also child ratings, indicated by c.

h2TabM <- 
  fread(paste0("Tables_GenSem/ChPsy/mvLDSC_ChPsyGenFscAsr_h2FullTable.txt"),header=T) %>% 
  distinct(Trait,.keep_all = T) %>% 
  mutate(Trait=str_replace_all(Trait,"Q0","_Q0")) %>% 
  mutate(Trait=str_replace_all(Trait,"Q1","_Q1")) %>% 
  separate(Trait,c("Trait","Age"),"_") %>% 
  mutate(alpha = ifelse(h2<0 | h2-SE<0, 1, 1)) %>% 
  mutate(Type = paste("GWA")) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  mutate(Age=str_replace_all(Age,"Q04","6mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q05","18mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q06","3yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q07","5yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q09","8yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>% 
  mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>% 
  filter(!Trait %in% excl)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h2TabV <- 
  fread(paste0("Tables_GenSem/ChPsy/mvLDSC_ChPsyGenQrsAsr_h2FullTable.txt"),header=T) %>%
  distinct(Trait,.keep_all = T) %>% 
  mutate(Trait=str_replace_all(Trait,"Q0","_Q0")) %>% 
  mutate(Trait=str_replace_all(Trait,"Q1","_Q1")) %>% 
  separate(Trait,c("Trait","Age")) %>% 
  mutate(alpha = ifelse(h2<0 | h2-SE<0, 1, 1))%>% 
  mutate(Type = paste("vGWA")) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  mutate(Age=str_replace_all(Age,"Q04","6mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q05","18mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q06","3yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q07","5yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q09","8yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>% 
  mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>% 
  filter(!Trait %in% excl)

h2Tab1 <- rbind(h2TabM,h2TabV)

h2Tab2 <- 
  h2TabM %>% 
    select(Trait,Age,h2_mGWA=h2,SE_mGWA=SE,Z_mGWA=Z, I_mGWA=Int, X2_mGWA=MeanChi2, L_mGWA=Lambda) %>% 
  full_join(
    h2TabV %>% 
      select(Trait,Age,h2_vGWA=h2,SE_vGWA=SE,Z_vGWA=Z, I_vGWA=Int, X2_vGWA=MeanChi2, L_vGWA=Lambda)
    )


N <- fread("Tables_GenSem/ChPsy/ChPsyGenQrsAsr_SampleSizes.txt",header=F) %>% 
  mutate(V1=str_replace_all(V1,"_RegenieMQt.Munged.sumstats.gz","")) %>% 
  mutate(V1=str_replace_all(V1,"_RegenieVQt.Munged.sumstats.gz","")) %>% 
  mutate(V1=str_replace_all(V1,"GwaSumStats/ChPsyGenQrsAsr/","")) %>% 
  mutate(V1=str_replace_all(V1,"Q", "_Q")) %>% 
  separate(V1, c("Trait","Age")) %>% 
  select(Trait, Age, N=V2) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  mutate(Age=str_replace_all(Age,"Q04","6mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q05","18mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q06","3yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q07","5yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q09","8yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>% 
  mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>% 
  filter(!Trait %in% excl)


h2_TabFull <- fread("Tables_Phe/xIRT_ModelStats.tsv") %>%
  mutate(Trait=paste0(Scale,Subscale)) %>% 
  select(Trait,Age=Wave,everything()) %>% 
  mutate(Rater=ifelse(grepl(x=Age,"C"),"Child","Mother")) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  mutate(Age=str_replace_all(Age,"Q04","6mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q05","18mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q06","3yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q07","5yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q09","8yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>% 
  mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>% 
  select(!matches("Scale"),!matches("Subscale"),!matches("_"),-AIC) %>% 
  inner_join(h2Tab2) %>% inner_join(N) %>% 
  filter(!Trait %in% excl) %>% 
  mutate(Age=factor(Age,levels=c("6mo","18mo","3yr","5yr","8yr","14yr(m)","14yr(c)"),
                        labels=c("6 m/o","18 m/o","3 y/o","5 y/o","8 y/o","14 y/o","14 y/o")
                         )) %>% 
  mutate(
    h2_mGWA=round(h2_mGWA,2),h2_vGWA=round(h2_vGWA,2),
    SE_mGWA=round(SE_mGWA,2),SE_vGWA=round(SE_vGWA,2),
    lSE_vGWA=round(h2_vGWA-SE_vGWA,2),uSE_vGWA=round(h2_vGWA+SE_vGWA,2),
    lSE_mGWA=round(h2_mGWA-SE_mGWA,2), uSE_mGWA=round(h2_mGWA+SE_mGWA,2)
         )

fwrite(h2_TabFull,"Tables_GenSem/ChPsy/ChPsyh2_FullTable.txt",col.names=T)

meanh2V  <- round(mean(h2_TabFull$h2_vGWA),2)
meanZV   <- round(mean(h2_TabFull$Z_vGWA),2)
meanh2Vz <- round(mean(h2_TabFull %>% filter(Z_vGWA>3) %>% .$h2_vGWA),2)
meanh2m  <- round(mean(h2_TabFull$h2_mGWA),2)
meanZm   <- round(mean(h2_TabFull$Z_mGWA),2)
meanh2mz <- round(mean(h2_TabFull %>% filter(Z_mGWA>3) %>% .$h2_mGWA),2)

Figure 3

Relationship between vGWA (y-axis) and mGWA (x-axis) SNP-based heritability estimated using LDSC from summary statistics following GWA and vGWA of childhood behavioural traits in the MoBa sample

Figure_3 <-
  ggplot(
    h2_TabFull,
    aes(x=h2_mGWA,y=h2_vGWA)) +
  geom_linerange(aes(ymin=round(h2_vGWA-SE_vGWA,2),ymax=round(h2_vGWA+SE_vGWA,2),color=Trait)) +
  geom_linerange(aes(xmin=round(h2_mGWA-SE_mGWA,2),xmax=round(h2_mGWA+SE_mGWA,2),color=Trait)) +
  scale_x_continuous(limits=c(-0.1,0.2),breaks=c(seq(-.1,.2,0.1)))+
  scale_y_continuous(limits=c(-0.1,0.2),breaks=c(seq(-.1,.2,0.1)))+
  geom_hline(yintercept=0,colour="limegreen",linetype="dotted") +
  geom_vline(xintercept=0,colour="limegreen",linetype="dotted") +
  facet_wrap(facets=vars(Age),ncol=3) + 
  theme_minimal() +
  scale_colour_viridis_d(option="plasma") +
  guides(color="none")

ggplotly(Figure_3, tooltip = c("Trait","h2_mGWA","SE_mGWA","h2_vGWA","SE_vGWA"))

Scatter plots display the point estimates (intersection) and standard errors for h²_mGWA (x-axis) h²_vGWA (y-axis) from all MoBa analyses. The plot is stratified by Age of MoBa participants. MoBa traits are labelled in the hover over boxes (interactive version) and can be viewed in S.Table 6.

S.Table 6

SNP-based heritability estimated using LDSC from summary statistics following mGWA and vGWA of childhood behavioural traits in the MoBa sample

datatable(
  h2_TabFull %>% select(Trait,Age,Rater,h2_mGWA,SE_mGWA,Z_mGWA,h2_vGWA,SE_vGWA,Z_vGWA),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

See table one for description of each trait. 2 Time-points named to reflect the average age of the children: 6mo, 6 months old; 18mo, 18 months old; 3yr, 3 years old; 5yr, 5 years old; 8yr, 8 years old; 14yr, 14 years old. Outcomes measured at all time-points are mother-rated, apart from at the 14yr time-point, where there are also child ratings, indicated by c.

Mean values:h2_vGWA = 0.02, Z = 1.04 h2_mGWA = 0.06, Z = 3.67 when Z>3:h2_vGWA = 0.06; h2_mGWA = * 0.08

datatable(
  h2_TabFull %>% select(Trait,Age,Rater,I_mGWA,X2_mGWA,L_mGWA,I_vGWA,X2_vGWA,L_vGWA),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

S.Figure 3

Relationship between SNP-based heritability (y-axis) estimated using LDSC from summary statistics following GWA and vGWA of childhood behavioural traits in the MoBa sample and psychometric scale reliability (x-axis)

Rel_h2 <-
  ggplot(
    h2_TabFull %>% 
      select(Trait, vGWA=h2_vGWA, mGWA=h2_mGWA, Reliability) %>% 
      pivot_longer(!c(Trait,Reliability),names_to="Effect",values_to="h2"),
    aes(x=Reliability, y=h2, group=Effect, colour=Effect,shape=Effect,label=Trait)) +
  geom_point() +
  geom_smooth(method=lm,se=FALSE) +
  scale_colour_viridis_d(option="viridis",direction=-1,begin=0.4,end=0.8)+
  theme_minimal() +
  directlabels::geom_dl(aes(label=Effect),method="smart.grid")+
  guides(color="none",shape="none")#+ 
  # ylab(bquote({h^2}))


plotly::ggplotly(Rel_h2, tooltip = c("Reliability","h2","Effect","Trait"))

S.Figure 3 shows the relationship between SNP-based heritability and scale reliability (Cronbach’s alpha) stratified by the type of analysis (GWA versus vGWA). Reliability, measured using Cronbach’s alpha, is the average correlation between each item and the total scale. For GWA, h2SNP increases with the reliability of the measure, whereas for vGWA, h2SNP decreases as the reliability of the measure increases.

S.Figure 4

SNP based genetic correlations estimated using LDSC from summary statistics following GWA (left) and vGWA (right) of childhood behavioural traits in the MoBa sample

A. mGWA

mCorMat <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenFscAsr_rG_est.csv") %>% as.matrix()
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Ch","")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q04","_6mo")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q05","_18mo")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q06","_3yr")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q07","_5yr")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q09","_8yr")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat), "Q11C","_14yr(c)")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat), "Q11M","_14yr(m)")
row.names(mCorMat) <- colnames(mCorMat)

mCorTab <- round(mCorMat,2)
# mCorTab[upper.tri(mCorTab)] <- NA
mCorTab <- reshape2::melt(mCorTab)
names(mCorTab)[3] <- "rg"

ErrMat  <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenFscAsr_rG_err.csv") %>% as.matrix()
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q04","_6mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q05","_18mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q06","_3yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q07","_5yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q09","_8yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11C","_14yr(c)")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11M","_14yr(m)")
row.names(ErrMat) <- colnames(ErrMat)

ErrTab  <- round(ErrMat,3)
# ErrTab[upper.tri(ErrTab)] <- NA
ErrTab  <- reshape2::melt(ErrTab)
names(ErrTab)[3] <- "SE"

mCorTab <- mCorTab %>% full_join(ErrTab,by=c("Var1","Var2"))

mCorTab$Z <- mCorTab$rg/mCorTab$SE

Pval <- fread("Tables_GenSem/ChPsy/mvLDSC_ChPsyGenFscAsr_rG_pval.txt")%>% select(Var1=V1,Var2=V2,Pval=V3) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieVQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieMQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),str_remove_all,"Ch")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q04","_6mo")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q05","_18mo")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q06","_3yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q07","_5yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q09","_8yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q11C","_14yr(c)")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q11M","_14yr(m)")) %>% 
  mutate(across(c(Var1,Var2),as.factor))
  
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")

mCorTab <- mCorTab %>% full_join(Pval,by=c("Var1","Var2")) %>% na.omit()  %>% 
  filter(!grepl(x=Var1,pattern=paste0(excl,collapse="|"))) %>% 
  filter(!grepl(x=Var2,pattern=paste0(excl,collapse="|"))) %>% 
  mutate(Z=round(Z,1))


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


mCor  <- mCorTab %>% select(Var1,Var2,rg,SE,Z)
vars  <- unique(c(levels(mCor$Var1),levels(mCor$Var2)))

excl <- c(
  "IcqFull","CbclFull","EasFull","CfqFull","ItseaFull","ScqFull","SdqFull",
  "CccFull","SprakFull","NhpicFull","RsdbdFull","PeqFull","SbFull",
  "SclFull","SppaFull","CbclAdhd"
)

vars  <- vars[!grepl(paste0(excl,collapse="|"),vars)]

combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos  <- data.frame(V1=vars,V2=vars) 
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")

mCor      <- 
  bind_rows(
    mCorTab %>% select(Var1,Var2,rg),
    mCorTab %>% select(Var1=Var2,Var2=Var1,rg)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(rg=ifelse(Var1==Var2,1,rg)) %>% 
  pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
  column_to_rownames("Var2") %>% 
  as.matrix()

mCor[is.na(mCor)] <- 0

mCor[mCor>1] <- 1
mCor[mCor<(-1)] <- (-1)

Err  <- 
  bind_rows(
    mCorTab %>% select(Var1,Var2,SE),
    mCorTab %>% select(Var1=Var2,Var2=Var1,SE)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(SE=ifelse(Var1==Var2,0,SE)) %>% 
 pivot_wider(names_from=Var1,values_from=SE) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Err[is.na(Err)] <- 0

Zsc  <- 
  bind_rows(
    mCorTab %>% select(Var1,Var2,Z),
    mCorTab %>% select(Var1=Var2,Var2=Var1,Z)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(Z=ifelse(Var1==Var2,2,Z)) %>% 
 pivot_wider(names_from=Var1,values_from=Z) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Zsc[is.na(Zsc)] <- 0
Zsc  <- abs(Zsc)


ChmCorplotly<-
    heatmaply::heatmaply_cor(
      mCor,
      node_type = "scatter",
      point_size_mat=Zsc,
      point_size_name="|Z|",
      label_names = c("x", "y", "rg"),
      fontsize_col=8,
      fontsize_row=8,
      grid_size = 0.001,
      dendrogram = "both",
      subplot_heights=c(0.06, 0.94),
      subplot_widths=c(0.94,0.06),
      hide_colorbar=T
    ) 

ChmCorplotly$height <- nrow(mCor)*14
ChmCorplotly$width  <- nrow(mCor)*14
ChmCorplotly$x$data[[4]]$marker$size <- (ChmCorplotly$x$data[[4]]$marker$size/2.5)

ChmCorplotly

Genetic correlation matrices for mGWA (left) and vGWA (right) shows the degree to which genetic effects are shared across the phenotypes. The colour intensity represents the direction and strength of the genetic correlation (red = positive correlation, white = no correlation, blue = negative correlation) and the size of each point represents the z-statistic (rg / se). For mGWA (left), the correlation structure is more dense than for vGWA (right). This is likely because of differences in statistical power for detecting mean and variance effects - as a much smaller proportion of vGWA summary statistics had statistically significant h2 estimates.

B. vGWA

vCorMat <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenQrsAsr_rG_est.csv") %>% as.matrix()
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"_RegenieVQt.Munged.sumstats.gz","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Ch","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q04","6mo")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q05","18mo")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q06","3yr")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q07","5yr")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q09","8yr")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q11C","14yr(c)")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q11M","14yr(m)")
row.names(vCorMat) <- colnames(vCorMat)

vCorTab <- round(vCorMat,2)
# vCorTab[upper.tri(vCorTab)] <- NA
vCorTab <- reshape2::melt(vCorTab)
names(vCorTab)[3] <- "rg"

ErrMat  <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenQrsAsr_rG_err.csv") %>% as.matrix()
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieVQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q04","6mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q05","18mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q06","3yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q07","5yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q09","8yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11C","14yr(c)")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11M","14yr(m)")
row.names(ErrMat) <- colnames(ErrMat)

ErrTab  <- round(ErrMat,3)
# ErrTab[upper.tri(ErrTab)] <- NA
ErrTab  <- reshape2::melt(ErrTab)
names(ErrTab)[3] <- "SE"

vCorTab <- vCorTab %>% full_join(ErrTab,by=c("Var1","Var2"))

vCorTab$Z <- vCorTab$rg/vCorTab$SE

Pval <- fread("Tables_GenSem/ChPsy/mvLDSC_ChPsyGenQrsAsr_rG_pval.txt")%>% 
  select(Var1=V1,Var2=V2,Pval=V3) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieVQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieMQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),as.factor)) %>% 
  mutate(across(c(Var1,Var2),str_remove_all,"Ch")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q04","6mo")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q05","18mo")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q06","3yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q07","5yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q09","8yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q11C","14yr(c)")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q11M","14yr(m)")) %>% 
  mutate(across(c(Var1,Var2),as.factor))

vCorTab <- vCorTab %>% full_join(Pval,by=c("Var1","Var2")) %>% na.omit()%>% 
  mutate(across(c(Var1,Var2),as.character))%>% 
  filter(!grepl(x=Var1,pattern=paste0(excl,collapse="|"))) %>% 
  filter(!grepl(x=Var2,pattern=paste0(excl,collapse="|"))) %>% 
  mutate(Z=round(Z,1))

## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

vCor  <- vCorTab %>% select(Var1,Var2,rg,SE,Z)
vars  <- unique(c(vCor$Var1,vCor$Var2))

excl <- c(
  "IcqFull","CbclFull","EasFull","CfqFull","ItseaFull","ScqFull","SdqFull",
  "CccFull","SprakFull","NhpicFull","RsdbdFull","PeqFull","SbFull",
  "SclFull","SppaFull","CbclAdhd"
)

vars  <- vars[!grepl(paste0(excl,collapse="|"),vars)]

combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos  <- data.frame(V1=vars,V2=vars) 
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")

vCor      <- 
  bind_rows(
    vCorTab %>% select(Var1,Var2,rg),
    vCorTab %>% select(Var1=Var2,Var2=Var1,rg)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(rg=ifelse(Var1==Var2,1,rg)) %>% 
  pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
  column_to_rownames("Var2") %>% 
  as.matrix()

vCor[is.na(vCor)] <- 0

vCor[vCor>1] <- 1
vCor[vCor<(-1)] <- (-1)

Err  <- 
  bind_rows(
    vCorTab %>% select(Var1,Var2,SE),
    vCorTab %>% select(Var1=Var2,Var2=Var1,SE)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(SE=ifelse(Var1==Var2,0,SE)) %>% 
 pivot_wider(names_from=Var1,values_from=SE) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Err[is.na(Err)] <- 0

Zsc  <- 
  bind_rows(
    vCorTab %>% select(Var1,Var2,Z),
    vCorTab %>% select(Var1=Var2,Var2=Var1,Z)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(Z=ifelse(Var1==Var2,2,Z)) %>% 
 pivot_wider(names_from=Var1,values_from=Z) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Zsc[is.na(Zsc)] <- 0
Zsc  <- abs(Zsc)

ChvCorplotly<-
    heatmaply::heatmaply_cor(
      as.data.frame(vCor),
      node_type = "scatter",
      # point_size_mat = -log(Pval),
      # point_size_name = "-logP",
      point_size_mat=Zsc,
      point_size_name="|Z|",
      label_names = c("x", "y", "rg"),
      fontsize_col=8,
      fontsize_row=8,
      grid_size = 0.001,
      dendrogram = "both",
      subplot_heights=c(0.06, 0.94),
      subplot_widths=c(0.94,0.06),
      hide_colorbar=T
    ) 

ChvCorplotly$height <- nrow(vCor)*17
ChvCorplotly$width  <- nrow(vCor)*17
ChvCorplotly$x$data[[4]]$marker$size <- (ChvCorplotly$x$data[[4]]$marker$size/2)

ChvCorplotly

Genetic correlation matrices for mGWA (left) and vGWA (right) shows the degree to which genetic effects are shared across the phenotypes. The colour intensity represents the direction and strength of the genetic correlation (red = positive correlation, white = no correlation, blue = negative correlation) and the size of each point represents the z-statistic (rg / se). For mGWA (left), the correlation structure is more dense than for vGWA (right). This is likely because of differences in statistical power for detecting mean and variance effects - as a much smaller proportion of vGWA summary statistics had statistically significant h2 estimates.

S.Table 7

SNP-based genetic correlations estimated using LDSC from summary statistics following GWA and vGWA of childhood behavioural traits in the MoBa sample

A. mGWA

RgM_N  <- nrow(mCorTab %>% filter(Z > 3))
meanMRgZ  <- round(mean(abs(mCorTab$Z)),2)

datatable(
  mCorTab %>% select(Var1,Var2,rg,SE,Pval,Z),
  filter='bottom',rownames=F,extensions='Buttons',
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

B. vGWA

meanVRgZ  <- round(mean(abs(vCorTab$Z)),2)
RgV_N  <- nrow(vCorTab %>% filter(Z > 3))

datatable(
  vCorTab %>% select(Var1,Var2,rg,SE,Pval,Z),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

Supplement 2: Analysis of height and weight

S.Table 8

Genome-wide significant SNPs identified in mGWA & vGWA analyses (independent loci r²<.1)

SnpTab <- fread("Tables_Snps/ChAnth/TopSnps_Clumped.tsv") %>%
  mutate(af = round(af,2)) %>% 
  mutate(mP = formatC(mP, format = "e", digits = 2)) %>% 
  mutate(vP = formatC(vP, format = "e", digits = 2)) %>%  
  mutate(across(c(Trait),str_remove_all,"Ch"))

DT::datatable(
  SnpTab,
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

CorP <- cor.test(as.numeric(SnpTab$mP),as.numeric(SnpTab$vP))

Correlation between p-values for mGWA and vGWA: Pearson’s r = 2.68; p = 7.50e-03

S.Figure 5

PheWas associations across all independent (r<.1) MoBa vGWA (top) and GWA (bottom) associations (p<1e-5)

vPheWasTab <- fread("Tables_Snps/ChAnth/Anth_vPheWasTab.tsv")
mPheWasTab <- fread("Tables_Snps/ChAnth/Anth_mPheWasTab.tsv")
PheWasTab3 <- vPheWasTab %>% bind_rows(mPheWasTab)

PheWasLab <- 
  PheWasTab3 %>% 
  arrange(PheWasP) %>% distinct(PheWasTrait,MoBaTrait,Age,.keep_all=T) %>% 
  select(PheWasTrait,MoBaTrait,Age,PheWasP) %>%
  mutate(Label=str_replace_all(PheWasTrait, "\\s*\\([^\\)]+\\)", "")) %>% 
  mutate(Label=str_replace_all(Label, "[[:punct:]]", "")) %>% 
  mutate(Label=str_to_title(Label)) %>% 
  mutate(Label=str_replace_all(Label, " ", "")) %>% 
  mutate(Label=str_replace_all(Label, "CellCount", "Count")) %>% 
  # distinct(Label,MoBaTrait,.keep_all=T) %>% 
  inner_join(PheWasTab3,c("MoBaTrait","PheWasTrait","PheWasP","Age")) %>% 
  select(chr,rsid,MoBaTrait,Age,PheWasTrait=Label,PheWasP,vP,mP) %>% distinct() %>% 
  mutate(MoBaTrait=str_replace_all(MoBaTrait, "C_", " ")) %>% 
  mutate(MoBaTrait=str_replace_all(MoBaTrait, "Full", "")) %>% 
  mutate(MoBaTrait=str_replace_all(MoBaTrait, "([a-z])([A-Z])", "\\1 \\2")) %>% 
  mutate(MoBaTrait=str_replace_all(MoBaTrait, "([a-z])([0-9])", "\\1 \\2"))%>% 
  mutate(Age = factor(Age, levels=c("18mo","03yo","05yo","07yo","08yo","14yo")))

PheWasLab <- 
  bind_rows(
  PheWasLab %>% filter(vP<1e-5) %>% rename("GwaP"=vP) %>% mutate(Effect="vQTL"),
  PheWasLab %>% filter(mP<1e-5) %>% rename("GwaP"=mP) %>% mutate(Effect="mQTL")
  ) %>% 
  mutate(GwaP=round(-log10(GwaP),1),PheWasP=round(-log10(PheWasP),1)) 


PheWasLab$mP <- formatC(PheWasLab$mP, format = "e", digits = 2)
PheWasLab$vP <- formatC(PheWasLab$vP, format = "e", digits = 2)


PheWasPlot <- 
  ggplot(PheWasLab,
    aes(y=PheWasP,x=GwaP,shape=Effect,colour=MoBaTrait,label=PheWasTrait))+ 
    geom_point()+theme_minimal()+
    facet_grid(cols=vars(Age),rows=vars(Effect),drop=T)+
    scale_colour_viridis_d(begin=.2,end=.8)+
    scale_shape_manual(values=c(0:14))+
    theme(legend.position="none")

ggplotly(PheWasPlot)

S.Table 9

Results from a phenome wide association scan of publicly available datasets

DT::datatable(
  PheWasLab %>% select(chr,rsid,MoBaTrait,Age,PheWasTrait,PheWasP,vP,mP),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

h2TabM <- 
  fread(paste0("Tables_GenSem/ChAnth/mvLDSC_ChAnth_h2FullTable.txt"),header=T) %>% 
  mutate(across(c(h2, SE, Z), ~round(.x,2)))  %>% 
  distinct(Trait,.keep_all = T) %>% 
  separate(Trait,c("Trait","Age","Adj"),"_") %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  filter(!(Trait=="Weight"&is.na(Adj))) %>% 
  mutate(Type = paste("mGWA")) %>% 
  select(-Adj)

h2TabV <- 
  fread(paste0("Tables_GenSem/ChAnth/mvLDSC_ChAnthQrs_h2FullTable.txt"),header=T) %>% 
  mutate(across(c(h2, SE, Z), ~round(.x,2)))  %>% 
  distinct(Trait,.keep_all = T) %>% 
  separate(Trait,c("Trait","Age","Adj"),"_") %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  filter(!(Trait=="Weight"&is.na(Adj))) %>% 
  mutate(Type = paste("vGWA")) %>% 
  select(-Adj)

h2Tab1 <- rbind(h2TabM,h2TabV)

h2Tab2 <- 
  h2TabM %>% 
  select(Trait,Age,h2_mGWA=h2,SE_mGWA=SE,Z_mGWA=Z) %>% 
  full_join(h2TabV %>% select(Trait,Age,h2_vGWA=h2,SE_vGWA=SE,Z_vGWA=Z))

N <- fread("./Tables_GenSem/ChAnth/ChAnthQrs_SampleSizes.txt",header=F) %>% 
  mutate(V1=str_replace_all(V1,"_RegenieMQt.Munged.sumstats.gz","")) %>% 
  mutate(V1=str_replace_all(V1,"_RegenieVQt.Munged.sumstats.gz","")) %>% 
  mutate(V1=str_replace_all(V1,"GwaSumStats/ChAnthQrs/","")) %>% 
  # mutate(V1=str_replace_all(V1,"([a-z])([A-Z])", "\\1 \\2")) %>% 
  separate(V1, c("Trait","Age","Adj")) %>% 
  select(Trait, Age, N=V2) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch",""))

meanh2V  <- round(mean(h2TabV$h2),2)
meanh2m  <- round(mean(h2TabM$h2),2)
meanZV   <- round(mean(h2TabV$Z),2)
meanZm   <- round(mean(h2TabM$Z),2)
meanh2Vz <- round(mean(h2TabV %>% filter(Z>3) %>% .$h2),2)
meanh2mz <- round(mean(h2TabM %>% filter(Z>3) %>% .$h2),2)

S.Table 10

SNP-based heritability estimated using LDSC from summary statistics following mGWA and vGWA of childhood behavioural traits in the MoBa sample

h2Tab2$Age <- factor(h2Tab2$Age,levels=c("18mo","03yo","05yo","07yo","08yo","14yo"))

datatable(
  h2Tab2, # %>% mutate(`|Zm|` = abs()),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

meanh2V  <- round(mean(h2TabV$h2),4)
maxh2V  <- round(max(h2TabV$h2),4)
meanh2m  <- round(mean(h2TabM$h2),4)
maxh2m  <- round(max(h2TabM$h2),4)

meanZV   <- round(mean(h2TabV$Z),2)
meanZm   <- round(mean(h2TabM$Z),2)
meanh2Vz <- round(mean(h2TabV %>% filter(Z>3) %>% .$h2),2)
meanh2mz <- round(mean(h2TabM %>% filter(Z>3) %>% .$h2),2)
h2V_N  <- nrow(h2TabV %>% filter(Z > 3))
h2M_N  <- nrow(h2TabM %>% filter(Z > 3))
meanXV   <- round(mean(h2TabV$MeanChi2),2)
meanXm   <- round(mean(h2TabM$MeanChi2),2)

Mean h2-vGWA = 0.0208, meanZ= 1.09; Mean h2-mGWA = 0.3008, meanZ= 10.11;
when Z>3: Mean h2-vGWA = NaN; Mean h2-mGWA = 0.32

S.Figure 6

Relationship between vGWA (y-axis) and mGWA (x-axis) SNP-based heritability estimated using LDSC from summary statistics following GWA and vGWA of childhood behavioural traits in the MoBa sample

Mh2_vGWAh2 <-
  ggplot(
    h2Tab2 %>% 
      mutate(
        h2_GWA=round(h2_mGWA,2),h2_vGWA=round(h2_vGWA,2),
        se_GWA=round(SE_mGWA,2),se_vGWA=round(SE_vGWA,2),
        lse_vGWA=round(h2_vGWA-SE_vGWA,2),use_vGWA=round(h2_vGWA+SE_vGWA,2),
        lse_GWA=round(h2_mGWA-SE_mGWA,2), use_GWA=round(h2_mGWA+SE_mGWA,2)
        ),
    aes(
      x=h2_GWA,y=h2_vGWA
      )) +
  geom_linerange(aes(ymin=round(h2_vGWA-SE_vGWA,2),ymax=round(h2_vGWA+SE_vGWA,2),color=Trait)) +
  geom_linerange(aes(xmin=round(h2_mGWA-SE_mGWA,2),xmax=round(h2_mGWA+SE_mGWA,2),color=Trait)) +
  scale_x_continuous(limits=c(-0.1,0.5),breaks=c(seq(0,0.5,0.1)))+
  scale_y_continuous(limits=c(-0.1,0.5),breaks=c(seq(0,0.5,0.1)))+
  geom_hline(yintercept=0,colour="limegreen",linetype="dotted") +
  geom_vline(xintercept=0,colour="limegreen",linetype="dotted") +
  # geom_smooth(method=lm,se=FALSE,colour="limegreen") +
  facet_wrap(facets=vars(Age),ncol=3) + #, scales = "free",space = "free")+
  theme_minimal() +
  scale_colour_viridis_d(option="plasma",begin=0.2,end=0.8)+
  guides(color="none")


ggplotly(Mh2_vGWAh2, tooltip = c("Trait","h2_GWA","se_GWA","h2_vGWA","se_vGWA"))

S.Figure 7

SNP based genetic correlations estimated using LDSC from summary statistics following GWA (left) and vGWA (right) of childhood behavioural traits in the MoBa sample

A. mGWA

mCorMat <- fread("Tables_GenSem/ChAnth/MvLDSC_ChAnth_rG_est.csv") %>% as.matrix()
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Ch","")
row.names(mCorMat) <- colnames(mCorMat)

KeepRows <- which(grepl(pattern="Height|Adj",rownames(mCorMat)))
KeepCols <- which(grepl(pattern="Height|Adj",colnames(mCorMat)))

mCorMat <- mCorMat[KeepRows,KeepCols]
colnames(mCorMat) <- str_remove_all(string=colnames(mCorMat),pattern="_hAdj")
row.names(mCorMat) <- colnames(mCorMat)

mCorTab <- round(mCorMat,2)
# mCorTab[upper.tri(mCorTab)] <- NA
mCorTab <- reshape2::melt(mCorTab)
names(mCorTab)[3] <- "rg"


ErrMat  <- fread("Tables_GenSem/ChAnth/MvLDSC_ChAnth_rG_err.csv") %>% as.matrix()
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
row.names(ErrMat) <- colnames(ErrMat)

ErrMat <- ErrMat[KeepRows,KeepCols]
colnames(ErrMat) <- str_remove_all(string=colnames(ErrMat),pattern="_hAdj")
row.names(ErrMat) <- colnames(ErrMat)

ErrTab  <- round(ErrMat,3)
# ErrTab[upper.tri(ErrTab)] <- NA
ErrTab  <- reshape2::melt(ErrTab)
names(ErrTab)[3] <- "SE"

mCorTab <- mCorTab %>% full_join(ErrTab,by=c("Var1","Var2"))

mCorTab$Z <- mCorTab$rg/mCorTab$SE

Pval <- fread("Tables_GenSem/ChAnth/mvLDSC_ChAnth_rG_pval.txt")%>% select(Var1=V1,Var2=V2,Pval=V3) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieVQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieMQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),str_remove_all,"Ch")) %>% 
  mutate(across(c(Var1,Var2),as.factor)) %>% 
  filter(grepl(pattern="Height|Adj",Var1))%>% 
  filter(grepl(pattern="Height|Adj",Var2)) %>% 
    mutate(across(c(Var1,Var2),str_remove_all,"_hAdj"))

  
mCorTab <- mCorTab %>% 
  full_join(Pval,by=c("Var1","Var2")) %>% na.omit()  %>% 
  mutate(Z=round(Z,1))


mCor  <- mCorTab %>% select(Var1,Var2,rg,SE,Z)

vars  <- unique(c(mCor$Var1,mCor$Var2))

combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos  <- data.frame(V1=vars,V2=vars) 
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")

mCor      <- 
  bind_rows(
    mCorTab %>% select(Var1,Var2,rg),
    mCorTab %>% select(Var1=Var2,Var2=Var1,rg)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(rg=ifelse(Var1==Var2,1,rg)) %>% 
  pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
  column_to_rownames("Var2") %>% 
  as.matrix()

mCor[is.na(mCor)] <- 0

mCor[mCor>1] <- 1
mCor[mCor<(-1)] <- (-1)

Err  <- 
  bind_rows(
    mCorTab %>% select(Var1,Var2,SE),
    mCorTab %>% select(Var1=Var2,Var2=Var1,SE)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(SE=ifelse(Var1==Var2,0,SE)) %>% 
 pivot_wider(names_from=Var1,values_from=SE) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Err[is.na(Err)] <- 0

Zsc  <- 
  bind_rows(
    mCorTab %>% select(Var1,Var2,Z),
    mCorTab %>% select(Var1=Var2,Var2=Var1,Z)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(Z=ifelse(Var1==Var2,2,Z)) %>% 
 pivot_wider(names_from=Var1,values_from=Z) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Zsc[is.na(Zsc)] <- 0
Zsc  <- abs(Zsc)


ChmCorplotly<-
    heatmaply::heatmaply_cor(
      mCor,
      node_type = "scatter",
      # point_size_mat = -log(Pval),
      # point_size_name = "-logP",
      point_size_mat=Zsc,
      point_size_name="|Z|",
      label_names = c("x", "y", "rg"),
      fontsize_col=8,
      fontsize_row=8,
      grid_size = 0.001,
      dendrogram = "both",
      subplot_heights=c(0.06, 0.94),
      subplot_widths=c(0.94,0.06),
      hide_colorbar=T
    ) 

# ChmCorplotly$height <- nrow(mCor)*15
# ChmCorplotly$width  <- nrow(mCor)*15
# ChmCorplotly$x$data[[4]]$marker$size <- (ChmCorplotly$x$data[[4]]$marker$size/2)

ChmCorplotly

B. vGWA

vCorMat <- fread("Tables_GenSem/ChAnth/MvLDSC_ChAnthQrs_rG_est.csv") %>% as.matrix()
colnames(vCorMat)  <- str_replace_all(colnames(vCorMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(vCorMat)  <- str_replace_all(colnames(vCorMat),"_RegenieVQt.Munged.sumstats.gz","")
colnames(vCorMat)  <- str_replace_all(colnames(vCorMat),"Ch","")
row.names(vCorMat) <- colnames(vCorMat)

KeepRows <- which(grepl(pattern="Height|Adj",rownames(vCorMat)))
KeepCols <- which(grepl(pattern="Height|Adj",colnames(vCorMat)))

vCorMat <- vCorMat[KeepRows,KeepCols]
colnames(vCorMat) <- str_remove_all(string=colnames(vCorMat),pattern="_hAdj")
row.names(vCorMat) <- colnames(vCorMat)

vCorTab <- round(vCorMat,2)
# vCorTab[upper.tri(vCorTab)] <- NA
vCorTab <- reshape2::melt(vCorTab)
names(vCorTab)[3] <- "rg"



ErrMat  <- fread("Tables_GenSem/ChAnth/MvLDSC_ChAnthQrs_rG_err.csv") %>% as.matrix()
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieVQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
row.names(ErrMat) <- colnames(ErrMat)

ErrMat <- ErrMat[KeepRows,KeepCols]
colnames(ErrMat) <- str_remove_all(string=colnames(ErrMat),pattern="_hAdj")
row.names(ErrMat) <- colnames(ErrMat)


ErrTab  <- round(ErrMat,3)
# ErrTab[upper.tri(ErrTab)] <- NA
ErrTab  <- reshape2::melt(ErrTab)
names(ErrTab)[3] <- "SE"

vCorTab <- vCorTab %>% full_join(ErrTab,by=c("Var1","Var2"))

vCorTab$Z <- vCorTab$rg/vCorTab$SE

Pval <- fread("Tables_GenSem/ChAnth/mvLDSC_ChAnthQrs_rG_pval.txt")%>% 
  select(Var1=V1,Var2=V2,Pval=V3) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieVQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieMQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),as.factor)) %>% 
  mutate(across(c(Var1,Var2),str_remove_all,"Ch")) %>% 
  mutate(across(c(Var1,Var2),as.factor)) %>% 
  filter(grepl(pattern="Height|Adj",Var1))%>% 
  filter(grepl(pattern="Height|Adj",Var2)) %>% 
  mutate(across(c(Var1,Var2),str_remove_all,"_hAdj"))



vCorTab <- vCorTab %>% full_join(Pval,by=c("Var1","Var2")) %>% na.omit()%>% 
  mutate(across(c(Var1,Var2),as.character))%>% 
  mutate(Z=round(Z,1))

vCor  <- vCorTab %>% select(Var1,Var2,rg,SE,Z)
vars  <- unique(c(vCor$Var1,vCor$Var2))

combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos  <- data.frame(V1=vars,V2=vars) 
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")

vCor      <- 
  bind_rows(
    vCorTab %>% select(Var1,Var2,rg),
    vCorTab %>% select(Var1=Var2,Var2=Var1,rg)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(rg=ifelse(Var1==Var2,1,rg)) %>% 
  pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
  column_to_rownames("Var2") %>% 
  as.matrix()

vCor[is.na(vCor)] <- 0

vCor[vCor>1] <- 1
vCor[vCor<(-1)] <- (-1)

Err  <- 
  bind_rows(
    vCorTab %>% select(Var1,Var2,SE),
    vCorTab %>% select(Var1=Var2,Var2=Var1,SE)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(SE=ifelse(Var1==Var2,0,SE)) %>% 
 pivot_wider(names_from=Var1,values_from=SE) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Err[is.na(Err)] <- 0

Zsc  <- 
  bind_rows(
    vCorTab %>% select(Var1,Var2,Z),
    vCorTab %>% select(Var1=Var2,Var2=Var1,Z)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(Z=ifelse(Var1==Var2,2,Z)) %>% 
 pivot_wider(names_from=Var1,values_from=Z) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Zsc[is.na(Zsc)] <- 0
Zsc  <- abs(Zsc)

ChvCorplotly<-
    heatmaply::heatmaply_cor(
      as.data.frame(vCor),
      node_type = "scatter",
      # point_size_mat = -log(Pval),
      # point_size_name = "-logP",
      point_size_mat=Zsc,
      point_size_name="|Z|",
      label_names = c("x", "y", "rg"),
      fontsize_col=8,
      fontsize_row=8,
      grid_size = 0.001,
      dendrogram = "both",
      subplot_heights=c(0.06, 0.94),
      subplot_widths=c(0.94,0.06),
      hide_colorbar=T
    ) 

ChvCorplotly

S.Table 11

SNP-based genetic correlations estimated using LDSC from summary statistics following GWA and vGWA of childhood behavioural traits in the MoBa sample

A. mGWA

datatable(
  mCorTab %>% select(Var1,Var2,rg,SE,Pval,Z) ,
  filter='bottom',rownames=F,extensions='Buttons',
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

B. vGWA

datatable(
  vCorTab %>% select(Var1,Var2,rg,SE,Pval,Z),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))